#pragma once
#include <vector>

// this class finds solution of the system of linear equations
// like this:
/*
	a11 * x1 + a12 * x2									        + a1n * xn + b1 = 0
	a12 * x1 + a22 * x2 + a23 * x3                                         + b2 = 0
	           a23 * x2 + a33 * x3 + a34 * x4                              + b3 = 0
	...
	          a(n-2)(n-1) * x(n-2) + a(n-1)(n-1) * x(n-1) + a(n-1)n * xn + b(n-1) = 0
	a1n * x1                                 + a(n-1)n * x(n-1) + ann * xn + bn = 0
-------------------------------------------------------------------------------------	
*/
// the matrix aij is symmetric, its elements are zero if abs(i-j) > 1 
// the diagonal vector is (a11, a22, a33, ... ann)
// the non-diagonal vector is (a12, a23, ... a(n-1)n, a1n)

// If the number of equations is 2, the system looks like:
/*
	a11 * x1 + (a12 + a23) * x2 + b1 = 0
	(a12 + a23) * x1 + a22 * x2 + b2 = 0
------------------------------------------------------------------------------------
*/

class SpecialEquationSolver
{
protected:
	typedef std::vector<float> FloatVector;

	FloatVector mDiagonal;
	FloatVector mNonDiagonal;
	FloatVector mB;
	FloatVector mX;
	int mEndIndex;
public:

	void Resize(int);

	void AddCoefficients(float diag, float nonDiag, float b);
	bool SolveLinearEquations();

	float GetX(int i) const;

	int GetIndexOfBiggestX() const;

protected:
	void RotateUp();
	void RotateDown();

private:

	// Recursive function; it finds the solution of the
	// reduced equations after some Gaussian elimination of variables.
	// Returns true on success.
	bool Step(int i, float diagFirst, float diagLast,
		float nonDiagLast, float bFirst, float bLast);
};

